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MAXIMUM RANGE THREE-DIMENSIONAL 
LIFTING PLANETARY ENTRY 

I. INTRODUCTION 


With the advent of the Space Shuttle, lifting reentry trajectories with 
medium lift-to-drag ratio (L/D) become of practical interest. Trajectories 
for this type of vehicle have, however, been investigated for a long time. 

The downrange capabilities were first discussed by Sanger [l] in 1933 and 
later on were studied by many investigators (see, e. g. , References 2 and 3). 
The trajectories are, in general, oscillatory and can be thought of as a 
quasisteady-glide (QSG) component (with y = 0) superimposed by a perturba- 
tion (with time-varying frequency and damping in the linear approximation). 

The QSG solution is an algebraic relationship between density, velocity, and lift 
coefficient, and this problem is readily solved for constant L/D. For small 
deviations from QSG, the equations may be linearized (Campbell [ 4]) or 
approximated by higher order expansions (Hanin[5]) to obtain frequency and 
damping of the actual trajectory. Recently Shi, et al. [6 and 7] using matched 
asymptotic expansions were able to achieve analytical approximate solutions 
to the nonlinear equations for a wide range of entry parameters. 

Cross-range investigations have appeared in the literature since about 
1960 [8 through 13]. Because of the increased dimension and nonlinearity 
of this problem, analytical solutions are even more sparse than for planar 
entry, and most investigations involve numerical computations. Eggers l 8] , 
Slye [9] , and Jackson tlO] consider a constant bank angle QSG maneuver. 
London [ 12] commented on the influence of the centrifugal term in the heading 
differential equation for these solutions. The equations obtained by Jackson 
for a spherical earth lend themselves to consider a bank angle program as a 
function of velocity and lateral range 111]. This minor circle turn is further 
discussed in References 13 through 16. 

In this report, a variational formulation for footprint maximization 
under QSG conditions is given. For a special case, some analytical approxi- 
mate integrals are obtained. For an entry vehicle with L/D = 1.4, QSG- 
footprint trajectories resulting from a numerical initial value search/trajec- 
tory optimization procedure are given. The perturbation equations around 
three-dimensional QSG trajectories are derived and discussed together with 
numerical solutions of non-QSG trajectories. Finally, the influence of a 
heating constraint on a maximum lateral range non-QSG reentry trajectory is 
discussed for a Space Shuttle orbiter-type vehicle. 



2. COORDINATE SYSTEM 


The equations of motion are written in a flight path oriented axis 
system, the x-axis of which is aligned with the velocity vector. The z-axis 
is downward in the vertical plane, and the origin of this right-handed system 
is based in the vehicle center of gravity (c. g. ) . The position of the vehicle 
c. g. with respect to the planet is specified by the downrange angle 0 , the 
cross-range angle A and the radial distance r = R+h (Fig. l). The 
heading angle relative to the initial conditions is designated by x , positive 
for right turns. 



Figure 1. Coordinate system. 


3. MATHEMATICAL MODEL AND SIMPLIFYING ASSUMPTIONS 

3. 1 Planet Model 

Because of the planet rotation, the oblateness, and the variation of 
the atmospheric properties with location and time, an actual entry trajectory 
depends on location, direction, and time of entry. For general entry trajec- 
tory investigations, a reduction of the possible parameter combinations is 
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often desirable. The dependency of the trajectory on entry time is eliminated 
by assuming 

(Ai) a stationary atmosphere, 

the dependency on entry point and direction by the following three assumptions: 
(A2) no planet rotation, 

(A3) no wind, 

(A4) planet and atmosphere are point symmetric (homogeneous on sphere- 
shells) . 

Assumption (A4) eliminates the oblateness effects and the latitudinal varia- 
tions of the atmosphere. For earth reentry with satellite velocity along the 
equator, assumption (A2) leads to errors of approximately 

6 percent in the aerodynamic velocity, 

12 percent in dynamic pressure, 

20 percent in the heating rate, 

initially. For meridional entry, the errors are much smaller. Best approxi- 
mations to actual trajectories can be expected with these assumptions for 
meridional entries in high latitudes. These assumptions about the planet 
considerably reduce the variational equations for reentry. 

Since the purpose of this report is to investigate basic properties 
of maximum range trajectories, they will be adopted along with others simpli- 
fying the vehicle model. The influence of oblateness and earth rotation on 
reentry trajectories is discussed in References 17, 18, and 19. The last 
two references also touch on the effect of assuming 

(A5) an exponential density law, p = P 0 exp(-^h) , 

which is convenient for analytical computations. 


3.2 Vehicle Model 

The vehicle is assumed to be 
(A 6) a point mass with instantaneous angular control. 
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The aerodynamic characteristics are assumed to be 
(A7) independent of Mach number and altitude 

and to obey the relationship (drag polar) 

(A8) C D =C D0+ kC L n . 

Assumption (A8) is a good representation for maximum range trajectories, 
which are flown in the vicinity of (L/D) . The consequence of assump- 

tion (A7) is to eliminate viscous interaction effects for high altitudes, which 
may have an appreciable influence on maximum range trajectories as shown 
in Reference 20. This can be accounted for partly by using a "medium" 
drag polar. Mach number effects become dominant in the supersonic and 
subsonic region which are not considered here. The trajectories are ended 
in 30 km altitude corresponding to about Ma = 3 . Although in the latter 
part of the trajectory, the dependency on the Mach number shows up, assump- 
tion (A7) is adopted throughout to reduce the number of parameters involved. 


3.3 Assumptions for Analytical Solution 

Further assumptions to obtain analytical approximate solutions are 

(A 9) small flight path angles: cos y = i , sin y = y 

(AiO) potential energy « kinetic energy: e « v/2 

(All) cylindrical planet: cos A = 1 

horizontal centrifugal component « aerodynamic turn force 

(A 12) Quasisteady-glide: y = 0 . 

The flight path angle will be determined from a differentiation of 
the QSG solution. 


3.4 Radiative Heating Constraint 

The heating constraint for a reentry vehicle is one of the driving 
factors in trajectory shaping. Medium- (like the proposed Space Shuttle 
orbiter) and high- L/D vehicles have a thermal protection system (TPS) 
consisting of reradiative and ablative elements at different parts of the vehicle 
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1 21 J . The higher the L/D, the more reradiative the TPS is going to be. If 
the time constants for the outer skin to heat up to the design limit is small 
compared to the rate at which the vehicle state (especially altitude, angle of 
attack, and velocity) changes, then a quasisteady approximation is valid. In 
this case, the heating constraint can be formulated as an algebraic equation 
<p(V, h,a) = 0, where V is the aerodynamic velocity, h is the altitude 
(representing air density) , and a is the angle of attack of the vehicle. It 
is assumed that the sideslip angle /3 is kept small so that its influence can 
be neglected. 


The data underlying the present representation are taken from 
Reference 22 for a limit temperature of 2000” F. Figure 2a shows a qualitative 
picture of typical altitude constraints due to dynamic pressure and kinetic 
heating. The effects of both velocity and lift coefficient are seen to be appre- 
ciable for the heating constraint. The heating constraint is more severe than 
the dynamic pressure boundary for hypersonic speeds (V > 2. 5 to 4. 5 km/s 
depending on the angle of attack). In the trajectory optimization program, on 
a boundary arc, the control C is required as a function of the state varia- 

bles altitude and velocity. The data are given in the form [22] 

tf>( h,V,a) = 0 . (3.1) 


Since the use of the lift coefficient C instead of angle of attack a reduces 

the computational workload to represent the aerodynamic characteristics of 
the vehicle, the relation C (a,h,V) is used to eliminate a . The trans- 

formed equation (3.1) is, then, represented in the form 


C LH = B i H i + iC LH ’ 1 = ^ 

where 

B i = g i i , J = 1.4 ; 

Hi = bh 2 /V 2 

H 2 = bh/V - H t 

H 3 = 1 - bh/V - H 2 

H 4 = V/(bh) - 2 + bh/V - H 3 

H 5 = V 2 /(bh) 2 - 3V/(bh) + 3 - bh/V - H 4 


(3.2) 


; (3.3) 
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and AC is an adjustment parameter; g 


LH 


1.3 


is a coefficient matrix 


0.110717 
-0. 672677 
0. 812241 
-3. 151267 
2.368095 


0. 834519 
2. 734170 
2.337815 
-13.621310 
19. 073400 


1.213679 
-0. 864369 
10.316280 
-40. 485500 
69. 869050 


-1.060833 
- 12. 100000 
22. 974860 
-57. 833330 
127. 777778 


(3.4) 


A 

h is a coordinate transformed altitude 
h = h/50tkm] - 1 


(3.5) 


to keep the coefficients closer to 1 . The constant b = 0. 095 serves the 
same purpose for the velocity. This form has been arrived at by intuition and 
trial. It is a bivariate polynomial of third order in the altitude and fourth 
order in velocity. With 


c, - 2 <1 - U * , , 

j=l ,J 

(3.6) 

the partial derivatives may be written 


ac LH /9V = b a/av(H ) 

(3.7) 

ac T „/ah = c. h. + b. a/ah(H.) 
LH' xi l ' l 

(3.8) 


From the form (3. l) one obtains for constant C the partial 

-LH 


ah 

av 



ac 


LH 


av 


ac 


LH 


dH 


The function generated in this way is shown in Figure 2b. 


(3.9) 
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The accuracy of the approximation is indicated in Figure 2b by dots 

representing the original input data into the curve fit procedure. These 

points had been obtained from Reference 22 through a crossplot. Except for 

the 70-km curve, the approximation is good considering the uncertainty in 

the original data. For maximum range and minimum energy loss trajectories 

which result in lift coefficients in the vicinity of maximum L/D (around 

C = 0. 2) , the 70-km curve will not be needed. Aside from that, the 
L 

deviation is on the safe side. 


4. EQUATIONS OF MOTION 


Under the assumptions made, 
D 


V = - g sin y - 


m 


X = - ~ cos y cos X tan A + 


y 


<? - i) 


cos y + 


mV 


_ V cos y 


9 = — 


r cos A 


cos x 


V 

A = ~ cos y sm x 
r 


r = R , = V sin y 
d 


the equations of motion may be written: 


L sin M 
mV cos y 


cos M 


U.1.1) 

(4. 1.2) 

(4.1.3) 

(4.1.4) 

(4.1.5) 

(4.1.6) 


where 



lift 


drag 


(4.1.7) 
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exponential atmosphere 

P = Po exp (-fS d h d ) 

gravity acceleration 

g = g 0 R 2 /(n + hj) 


(4.1.8) 


(4. 1. 9) 


The lift coefficient C and the bank angle M are controls. For analytical 

L 

solutions, a nondimensionalized form of these equations is advantageous. 
With the nondimensional time, 


and 


T = t\lg/l 


v = V 2 /(rg) , h = h /r 


(4.1.10) 


(4.1.11) 


these equations (4. 1. 1) to (4. 1. 7) become the following, where 


()=■—() 

dT 


a = Mr 

2m 


(4.1.12) 
(4. 1. 13) 


v = - \ fv~ 2 jja. exp ( — /3h) v + sin -yj 

X = \/v" j^a exp (-/3h) sin/ycos y - cos y cos x tan aJ (4. 1. 14) 

y = \l~v ^a exp ( -Ph) C ^ cos M - (l-v)/vcosyJ 

9 = \Tv cos y cos x/cos A 
A = n/v" cos y sin x 
h = sfv sin y 


(4.1.15) 

(4.1.16) 

(4.1.17) 

(4.1.18) 


The reciprocal atmospheric scale height now has approximately the value 
(for earth) 


P «* 930 


(4.1.19) 
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The solution to these equations for ma x imum range trajectories will 
now be approximated by superposition of two motion components: 

1. A QSG motion with equation (4. 1. 15) satisfied by y =0 

2. A small perturbation around the QSG trajectory. 


4.1 Qua si steady -glide Maximum Range 

4.1.1 VARIATIONAL QSG-EQUATIONS 


By introducing the specific energy 

e = v/2 + h (4. 1. 20) 


and adding equations (4. 1. 13) and (4. 1. 18) , the nondimensional altitude h 
becomes control-like, and, with v= 2(e - h) from equation (4.1.20) and 
applying assumption (A2) , the following set of equations are obtained for 
the analytical QSG solution: 


e = - <\/2(e - h) a exp (-/3h) 
X = s/2(e - h) £a exp (-/3h) 


0 = \l 2 (e - h) cos y cos x/cos A 
A = \/2(e - h) cos y sin x 
N/2(e - h) ^a exp (-/3h) cos M 


2(e - h) (4.1.21) 

sin M/ cos y - cos y cos x tan aJ 

(4.1.22) 

(4.1.23) 

(4.1.24) 

1 j- cos = 0 

(4.1.25) 


- | l/[2(e - h)] - 


These are four differential equations and one algebraic one, equation (4. 1.25) , 
which, however, has to be satisfied all along the trajectory. 

All initial conditions of the differential equations are considered to be 
given, while the initial altitude has to be found such that the payoff quantity 
obtains its extremum. This is an optimization problem with mixed side 
constraints (23J. The variational Hamiltonian is 
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H = ^ 2(e - h) < a exp ( -fih) -AC 


D 2(e - hj + C L ( 


“ — sin n + v cos M 
cos y 


) 


(- \ a + v° os A ) oos x + a a sin x - ” (2(6-10 ■ ‘) 


cos y) 


(4.1.26) 

and, since final time is open, the bracket { } is zero, v is a time-varying 
multiplier associated with the equality constraint (4. 1.25) and will be 
determined as a linear combination of the multipliers A [23] . 


For a moment, eliminate equation (4. 1.25) by introducing it into 
equations (4. 1.21) and (4. 1.22) in order to obtain a convenient expression 
for the optimal lift coefficient. Later on, however, one continues with the 
mixed side constraint formulation. One obtains 


e = 


- n/T 


1 - v 
E cos M 


cos y 


(4.1.27) 


where 


E = C l/ C D 

X = n/v~ ( tan — - cos y cos X tan a\ , (4.1.28) 

\cosy ) 

and the control-dependent part of the Hamiltonian in this formulation is with 
cos y « 1: 

H = \/v*<-A J 1 ~ ^ + A tan ju + . . .1 (4.1.29) 

) e E cos MX ) 

=1= 

3H /8C = 0 yields dE/dC =0 , i.e. , the lift-to-drag ratio 

L L 

(4.1.30) 

E has to be a maximum for a maximum range QSG trajectory. Since the 
aerodynamic drag polar is considered to be constant in the hypersonic range 
considered, E is a constant as well as C , the corresponding lift 

J-J-Lj 

coefficient. This information will be used when one returns, now, to the 
mixed formulation (4. 1. 26) . 
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controls C T , M and altitude h yields for maximum range ^ m ^ n : 


Looking for the extremum of the Hamiltonian with respect to the 
1 

'L 

Bank angle : 


tan V 


opt v cos y 


(4.1.31) 


or 


sin p = A /(-wcosy) , cos A* = iv/(-w) 
X 


w 


= ^v 2 + X^/cos 2 yj 


y 2 


Optimal lift coefficient : 
8C, 


D 


0C, 


-w 1 


A V E 

e 


'LE 


(4.1.32) 


(4.1.33) 


From equation (4.1.30) one concludes that w/(vA ) is a constant. Its 

e 

value is the inverse of the lift-to-drag ratio. 3H/3h = 0 furnishes the 
relationship 


a exp(-/3h) 


A e °DE 


Wv +2) + P w C 


LE 


2v 

~ ~r cos y - o 


(4.1.34) 


[a e X p(-») C D ] = v 2 [x e (/3v 2 +2) + /3Ew] • 

which with equation (4.1.33) reduces to 

fa exp(-/3h) C 1 = — ? • (4.1.35) 

L ° J opt A e V 
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Combining this with the QSG equation (4. 1.25) yields 


A 

v = -§■ v \l ( 1 - v) 
E 


(4.1.36) 


From equation (4. 1. 26) , the differential equations for the Lagrangian 
multipliers become (final time open, H = 0 ) 

X e = -«/2(e- h) |^a expl-Ph) C^-x^ + >- 2(e ‘ h)2 

and with equation (4. 1.35) 


A = 0 or A = const, 
e e 


A^ = - \f2(e - h) ^A@ - A^ sin A^ (-sinx) + A^, 


cos x 


0 


0 or A = const. 
0 


(4.1.37) 

(4. 1.38) 
(4. 1. 39) 


X. = 


-A — XT ( x ~ A_ sin a' 

A cos A \ X 0 


(4. 1.40) 


There are only two differential equations for the adjoints, and because of 
A g = const, the multiplier v is a function of v alone. 

The boundary conditions for maximum-range-QSG with a given kinetic 
energy decrease ^ and open final heading angle are 



att o 


at t f 

multiplier at t^ 

V 0 

= given 

= given 

^e„ = unknown = A 
f e 

Xo 

= 0 

> 


X f = open 

A =0 
X f 

©0 

= 0 


<fi = A + p9 = max 

CD 

II 

1 

*a 

Ao 

= 0 


II 

0 

(see section 4. 3. 3) 


A = -1 
A f 


h 0 , Mo , To have to be found such as to maximize the payoff quantity <p 
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From x = 0 and equation (4. 1.31) the equation for the final bank angle 
X i 

follows: 

M f = 0 . (4.1.41) 


Differentiating the QSG equation (4. 1,25) with respect to time and inserting 
equations (4. 1. 13) and (4. 1. 18) , an expression for the flight path angle is 
obtained: 


o 

_ 1 - v 2/E + \[y M sin M 

sin y cos ^ 2 + ( 1 - v) Pv 


(4.1.42) 


With equation (4.1.41) for the final time, this reduces to 
_ . -2(1 - v) 

Esmy i ’ 2 + 'llWW • 


(4.1.42a> 


Compared with results from numerical optimizations, this equation predicts 
the final flight path angle with good accuracy. It is the same as for planar 
entry. Except for the final time, the flight path angle is a function of the 
bank angle and its rate of change, however, 

4.1,2 SOME ANALYTICAL APPROXIMATE RESULTS FOR MAXIMUM 
CROSS-RANGE 

With X - 0 , which characterizes this case, and assumptions (A9) 
© 

to (All) , the variational QSG equations reduce to 
$ = -2 N/T^a exp(-^h) C^v 

X = n/’v a exp(-)Sh) C^ sin A* 

© = \fv cos x 

A = n/T sin x 


9 

X 


= 0 


0 

X 


- \Tv X. COS X 
A 


x e = 0 


(4.1.43) 


V. = 'sTv X cos X 
A X 
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(QSG) 


a exp(-/3h) C T v cos = (1 - v) /v 

with v = A vn/1 - v /E 
e 

where E = (L/D) = -vA /(\ 2 + v 2 V 2 

max e V X / 

Combining the last two lines with three equations for the four quantities p , 

A , v , and the product (A v) yields the surprising result 
X 6 

tan ix = n/v/( 1 - v) or sin ^ = n/v* . (4.1.44) 


J-lH. 


tan u , = A /v 
opt X 


C = C 

L opt LE 


Note that this simple relationship for the bank angle as a function of the 
velocity (m = arc sin of the ratio velocity to local satellite velocity) has 

T 

been obtained from the three relations 8H/3u = 0 , u = (C , n , h) , 

L 


together with the QSG condition without taking into account any boundary 
conditions. Therefore, this solution will be called "basic" and will be regarded 
with due suspicion. It may further be noted that this bank angle program has a 
square root relationship to the minor circle turn bank angle program [ il, 13 
through 16 ] . This solution is discussed below. 


4. 1.2. 1 "Basic" Solution . Differentiation of equation (4.1.44) with 
respect to time yields, with equations (4. 1. 13) and (4. 1.25) , 




-1 

E 



E sin y ^ 

\l 1 - v ) 


(4.1.45) 


Combining these results with equation (4. 1.42) , obtained by differentiating 
the QSG condition, leads to 


2 l_ (2 - v) \(1 - v 

E 2 + v[/3( 1 - v) - 1] 


and 


P = 


zi 

E 




2 - v ) 

2 + v[/3(l -v) - i]j 


(4.1.46) 


(4.1.45a) 


These results are plotted in Figure 3. It is seen that the flight path angle is 
small except for very small velocities. The bank angle rate is close to 
- i/E except for v close to 0 and 1 . 
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From the QSG condition, altitude as a function of velocity and lifting 
factor is 


h = 




(4.1.47) 


Inserting equations (4. i. 25), (4. 1.44) , and (4.1.46) into equation (4. 1. 13) 
yields 

* = f ' yv(1 " v) I 1 - 2+ vim"- V “7]} • <4 - 1 -‘ 

Neglecting the second term inside the brace [assumption (All), compare 
equation (4.1.45a) and Figure 3], this equation may be integrated to yield 
the final time: 




(4.1.49) 


Combining the first two differential equations of equation (4.1.43) yields, 
with equation (4. 1.44) , 



This equation neglects the centrifugal term in the heading differential equation 
so that for larger cross-ranges the actual heading angle will be smaller. 
Figure 4 shows how the inertial term grows with increasing lift-to-drag ratio. 
For L/D = 3 , the inertial summand averages about 15 percent of the aero- 
dynamic one. 

For orbital reentry with v 0 = 1 and v^ = 0.01 , equation (4.1.50) 

predicts a final heading of > ff/2 for vehicles with E > n/i. 8 = 1. 75 . 

Some results for other initial and final conditions and the vehicles parameters 
given in Table 1 are shown in column 8 of Table 2 . 
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Figure 4. Numerical check on assumption (Alla) for 
maximum lateral range trajectories. 


TABLE 1: VEHICLE CHARACTERISTICS 


E = L/D 

C D0 

k 

n 

C LE S / 

m 2 


kg_ 

1.4 

0.109 

1.45 

2.18 

0.001885 

2.22 

0.04 

1 . 

1.86 

0.000768 

2.38 

0.04 

0.95 

1.91 

0.0008 

3. 

0.028 

1.46 

2.21 

0.000612 


Wing Loading 
m/S /[kg/m 2 ] 

150 

250 

250 

250 


18 



TABLE 2. MANEUVER TIME AND FINAL HEADING 
FROM ANALYTICAL SOLUTIONS 


i 

2 

3 

* ' 

5 

6 

7 

8 





Maneuver times (t/sec) 







Cross-range 


Wo 

v o 

for n = 0 
at = 30km 

Downrange 
eq. (4.1.54) 

eq. (4.1.53) 

eq. (4.1.49) 

Numerical 

QSG 

eq. (4.1.50) 

1.40 

0.90850 

0.00822 

2024. 

1371. 

1338. 

1340. 

69. 

1.40 

0.95000 

0.00822 

2382. 

1471. 

1431. 

1458. 

71. 

1.40 

0.98180 

0.00822 

2968. 

1583. 

1534. 


72. 

2.22 

0.90850 

0.01990 

3118. 

2112. 

2030. 


103. 

2.22 

0.95000 

0.01990 

3684. 

2275. 

2178. 

2193. 

106. 

2.22 

0.98180 

0.01990 

4613. 

2460. 

2341. 

2382. 

108. 

3.00 

0.90850 

0.02490 

4171. 

2826. 

2702. 

2747. 

137. 

3.00 

0.95000 

0.02490 

4937. 

3049. 

2902. 

2880. 

140. 

3.00 

0.98180 

0.02490 

6192 

3302. 

3122. 


143. 


By assuming a linear bank angle program with time as suggested by 
equation (4.1. 45a) and Figure 3, the QSG differential equation for the velocity 
may be written as 

o _ _ 2 \fy ( 1 - y) 

v E cos (/u 0 - at) 

or 


dv 

n/7 (1 - v) 


2 dT 

E cos /u 0 cos ar + sin n 0 sin ar 


which integrates to 


In 


1 + n/v 


•1 f 


1 - n/v 



(4.1.51) 


. (4.1.52) 
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Requiring ^ = 0 at the final point leads to 


T = T 

f, L f, S 


h / In otn (j - f^J , 


(4. 1.53) 


where t„ c is Sanger's downrange maneuver time [1]: 
1> s 


r = 0.5E In 
i, b 


( 1 + n/vq) ( 1 - n/v) 
( 1 - n/v 0 ) ( 1 + n/v) 


(4.1.54) 


Columns 4 through 8 of Table 2 show maneuver times for maximum range 
trajectories ending at h^ = 30 km altitude. For conversion to dimensional 

time, an average value of Vr/I = 815. sec has been used. Sanger's down- 
range result [ 1 ] is given in column 4 . The next column contains the linear 
bank -angle -law cross -range maneuver time with Po taken from equation 
(4.1.44) using equation (4.1.53). The basic solution (4.1.49) yields column 
6, which is in surprising agreement with maneuver times obtained from opti- 
mizing numerically with the full set of nonlinear differential equations [equa- 
tions (4.1.1) to (4.1.9) and (4.3.1) through (4.3.9)] and QSG initial condi- 
tion h 0 and y 0 determined from equations (4.1.47) and (4.1.46), column 7. 
Except for the last case shown, the numerical results lie between those of 
columns 5 and 6. Maximum lateral range trajectories do have a final heading 

7 r 

angle in the vicinity of — relative to the entry direction. From this observa- 

tion and the numbers shown in column 8 for the final heading angle of the basic 
solution, the deviation of actual bank angle programs from the basic ones can 
be guessed. As will be discussed later in connection with the numerical 
results, optimal QSG maximum lateral range trajectories follow the basic 
bank angle program initially very closely and deviate from it toward the end, 
depending on L/D and the velocity range covered. 


A full analytical solution to the boundary value problem, equation 
(4. 1.43) , has not been found. Two integrals which check favorably with 
numerical results are the following: 


Combining the 

yields 


A and V 
X A 


differential equations in equation (4.1.43) 


A dA = -A a dA A 
XX A A 
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with the integral 


A 2 + A 2 = 1 , (4.1.55) 

X A 

here the integration constant has been determined from A = -1 and 

A f 

A = 0 . Despite the approximation made, this result agrees very well with 
X f 

-vierical solutions to the complete equations (4. 1. 1) to (4. 1. 6) as may be 
seen from Figure 5. 



Figure 5. Comparison of analytical approximate result 

A 2 + A 2 = 1 with numerical solutions. 

X A 
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o o 

Combining the © and A 


X 


equation using equation (4. i. 53) yields 



\J i - X 1 

X 


with the solution 

© = 0 O + arc sin A - arc sin A . (4.1.56) 

Xo X 

For ©o = 0 , the downrange at t^ , then, is 

©. = arc sin A . (4. 1. 56a) 

f Xo 

4.1.3 NUMERICAL QS G SOLU TIONS 

Even though oscillatory trajectories yield maximum range for non- 
equilibrium initial conditions, QSG trajectories are of special interest since 
they yield more uniform load distributions and passenger comfort without 
sacrificing much range. They need not, of course, be flown at maximum 
L/D if the mission does not require it. In this section, some numerically 
obtained QSG-footprint trajectories are described. 

To get an optimal three-dimensional QSG trajectory, initial conditions 
have to be matched with the unknown optimal control. This was done in an 
iteration loop on the initial conditions altitude and flight path angle around the 
optimization loop to determine the control. The optimization scheme used, 
a hybrid between the gradient and indirect methods, is especially suited for 
maximum range footprint computations [ 24] . The equations of motion 
underlying these computations are (4. 1. 1) to (4. 1. 9) with a drag polar of 
the form (A 8) . The adjoint equations are given in section 4. 3. 

Because of the computational workload involved, trajectories were 
considered to be reasonably converged QSG when, in the altitude plots, 
oscillations were hardly recognizable. Some of the trajectories could still 
be improved with more iterations. A good indication of optimal QSG is that 

the lift coefficient stays at C all the time. It is especially sensitive 

_LJd 

to non-QSG initial conditions. 
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Figure 6 shows cross-range versus downrange for maximum cross- 
range trajectories and three different vehicles (Table l) from different 
initial conditions (Table 3). 



0/deg 

Figure 6. Maximum lateral range QSG trajectories. 

TABLE 3. BOUNDARY CONDITIONS FOR MAXIMUM 
CROSS-RANGE QSG TRAJECTORIES 


L/D 

v 0 

To/0 

h 0 /km 

m °qsg/° 

Remarks 

1.4 

0.9085 

-0. 1729 

70.74 

72.4 

h^ = 30 km 


0. 95 

-0.2125 

73. 04 

77. 1 

QSG final conditions 

2.22 

0.95 

-0. 134 

66. 85 


see section 4. 3. 3 


0. 9818 

-0. 182 

70.37 

82.4 


3. 

0. 9085 

-0.0804 

63.59 




0. 95 

-0.0994 

65.28 
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A comparison of the (L/D = 2.22, v 0 = 0. 9818) - QSG trajectory 
with an oscillatory maximum cross-range trajectory with a flight path angle 
of y 0 = -1 deg in ho = 64 km altitude is given in Figure 7a. It is seen that 
the downrange is decreased by about 20 percent and that the cross-range 
remains the same. The corresponding controls are shown in Figures 7b and 
c. From Figure 7d, it may be seen that the oscillatory trajectory in an 
altitude-velocity plot centers around the QSG solution (see also section 4.2 on 
perturbations). Below V = 2 km/s , both trajectories are very similar. 




b. Bank angle time histories. 
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Figure 7. Maximum lateral range orbiter reentry 
L/D = 2.22 , m/F = 250 [kg/m 2 ]. 





c. Lift coefficient. 



V/(km/s) 


d. Altitude versus velocity 
Figure 7. (Concluded). 



The bank angle time histories in Figure 8 for v 0 = 0. 95 indicate the 
validity of the basic solution bank angle rate. Initially, it is very closely 
l/E , and it increases or decreases toward the end, depending on the lift-to- 
drag ratio E . For maximum cross-range, the final heading angle x^ is in 

the vicinity of 90 deg relative to entry. The analytical basic solution yields 
smaller values of x^ , equation (4. 1. 50) , for E = 1.4 and larger ones for 

E = 2.22 or 3 . Figure 9 shows how the adjustment in the bank angle 
'program is made by increasing the bank angle M at lower velocities for 
E = 1.4 and decreasing it for the higher values. If the initial velocity for 
E = 1.4 is lowered to 0. 9085 (dashed curve) , the bank angle is increased 
over the entire trajectory. 



Figure 8. Bank angle time histories for maximum lateral range QSG. 
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Figure 9. Comparison of bank angle programs between 
basic analytical approximate QSG solution and 
numerical solutions for maximum lateral 
range from same initial conditions. 

In Figure 10, QSG trajectories for different points on a footprint are 
given. It is seen that the bank angle program in the general case is not 
nearly linear with time. For maximum downrange, M is of course identically 
zero. 
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4.2 Perturbation Equations Around Quasisteady Glide 


For slightly perturbed initial conditions away from the QSG combina- 
tion, a good approximation to the actual trajectory may be obtained from 
.eries-developing the equations of motion around the QSG state as has been 
shown for the two-dimensional case in References 4 and 5. In this section, 
a similar approach is taken for the three-dimensional case based on the basic 
solution developed in section 4. 1. 2. 

With the state and control variable approximation 

x = x„ + <5x , u = u + 6u (4.2.1) 


where the index Q indicates the QSG solution, a Taylor series expansion of 
the differential equations (4. 1. 13) to (4. 1. 18) yields, with a truncation after 
the linear term in the perturbation, 


+ 0 = f ( x Q ’ u q) 


jif 

9x 


<5x + 


Q 


_8f 

au 


6u + 


Q 


(4.2.2) 


Removing the quasisteady glide component leads to the linear set of differen- 
tial equations with time-varying coefficients 

x = Ax + Bu , (4.2.3) 


where the 6-sign has been dropped for convenience and the coefficient 
matrices A and B are defined below [for x T = ( <5v , 6y , 6y , 6© , 6a) , 

u T = ^6 C l , , with 

b = a exp(-/3h) (4.2.4) 


and the QSG condition (4. 1.25) introduced into the matrix A , 
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0 Tf UJS q^o- 

0 rf soo q 


rl nis q 


0 

3 /Aq Z' 


(4.2.5) 


0 

0 

X SOO Ap 

0 

X tns A Ap - 

( a pz)/X *T S 

v SOD 

0 

V soo 


1121 a Ap - 

V SOO Apz 

V mrj X soo Ays. 

X tns Ap- 


X soo — I 

X soo 

V z soo 
X SOO AJ\- 

0 

V X UTS AJ\ 

Ap/ri ire; (A - i)gf- 

soo + ri rnq a ^ A Ap 

W/X 

0 

0 

0 

0 

_a p 

Ap/^ 9‘0 

0 

0 

0 

A/S./(A - T)tf- 

Ap/(A - T) A 





rl soo a 

a pz- 

L + x s °°a'\ : 

0 

0 

0 

A_p(A ~ 

\ ( A - T ) £ y 


v (4.2.6) 


4. 2. 1 FROZEN STATE APPROXIMATION 


Under the assumption that the variables on the reference trajectory 
are constant, momentary natural frequencies and damping of the perturbation 
motion can be computed. From the matrix A , equation (4.2. 5) , it is seen 
that the vertical motion is uncoupled from the rest (upper left quadrant) . It 
has the characteristic equation ( Laplace transform) : 


s 3 

(1 - v) 


+ s2 


(e^ - t ) + s 


JL 

E cos M 


- i 


)' 


E cos M (1 - v) 


+ /?v { y + [ (1 - v)y 2 - (i - 3v)]/(E cos m} = 0 


which over a wide region of v with y « i and 1 « fiv can be approximated 
by 


S 3 3 / 1 _ 3v \ 

TTTrt + sJ - sl3v + ( r - 0 • (4 - 2 - 7a) 

In Figure 11 the eigenvalues (period, damping) for the basic solution with 

cos M = n/T - v are plotted against the velocity parameter v using equation 
(4.2.7). 


30 




— ► V QSG 


a. Period and damping ratio of oscillatory component. 



V 


b. Unstable real root and oscillatory damping term. 
Figure li. Frozen state perturbation around basic QSG solution 


From Figure iia, it may be seen that the oscillation period, given in 
nondimensional time, is almost independent of the lift-to-drag ratio and 
decreases monotonically with deceleration. The damping ratio decreases 
with E = L/D increasing; as function of the velocity parameter v , it has 
a minimum around v = 0. 6 and assumes values around £ = 0. 4 for small v . 
Figure lib shows the unstable real root A* and the damping term of the 
oscillatory component. Though opposite in sign, they; follow the same trend 
characterized by a sharp increase in absolute value for v ^ 1/3 . For 
vehicles with higher L/D the subsidence mode is less unstable. 

4.2.2 NUMERICAL INT EGRATION OF PERTURBATION EQUATION 

For a linear bank angle law and an initial velocity corresponding to 
v 0 = 0.95 , three sets of equations have been integrated to check the approxi- 
mation of nonlinear oscillatory trajectories by a superposition of a QSG 
trajectory and linearized perturbation trajectories. The QSG trajectory was 
computed numerically using the same equations (4. 1. 13) to (4. 1. 18) , as for 
the non-QSG trajectory but initial conditions determined from the basic solu- 
tion. The full set of perturbation equations including the lateral motion but 
no control perturbation have been used. 

Figure 12 shows the altitude perturbation for a vehicle with E = 1.4 
and 6hfl = ±1.46 km and 6y 0 = -yQg G = 0.2125 deg. The oscillations 

for unchanged control do not center around <5h = 0 but around a displaced 
divergent curve when plotted over time (Fig. 12a). The period oscillation is 
well approximated by a medium value of the frozen state results from the 
previous section. Toward the end, the strong divergence predicted by the 
frozen state approximation becomes obvious. When the altitude is plotted 
over the velocity parameter v , however, all perturbed trajectories converge 
in the final part toward the QSG trajectory (Fig. 12b) . As seen from Figure 
7b, the optimized oscillatory lateral range trajectory oscillates around the 
basic solution. This shows that the control change due to non-QSG initial 
conditions is determined such that (see Figures 7b and c) the altitude oscil- 
lation takes place around the optimal QSG values. The lower curve in Figure 
12a (dashed) shows the difference between nonlinear non-QSG and QSG trajec- 
tory, and it demonstrates the good degree of approximation of the nonlinear 
solution achieved by the perturbation equations superimposed on a QSG trajec- 
tory. For larger initial perturbations and longer maneuver times, the quality 
of the approximation deteriorates. Figure 13 shows for 6y 0 = -0.2875 deg and 
E = 2.22 the discrepancy between linearized perturbation and nonlinear 
solution, which to a large extent stems from the first dip. In the nonlinear 
model, the first maximum in Ah is 20 percent lower than for the linearized 
equations, resulting in a time shift. The period of oscillation is slightly 
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I 



Figure 13. Comparison of perturbation solution to nonlinear 
results, E = 2. 22 , v 0 = 0. 95 , dy 0 = -0. 2875 deg , <5h 0 = 0 , 

smaller in the nonlinear case. The general tendency of the nonlinear solution 
is reflected correctly. The maximum deviation of dv occurs in the vicinity 
of v = 1/3 . 

The s im ple linear approximation used here seems to I.e valid for 
perturbations of a few tenths of a degree in flight path angle and a few kilo- 
meters in altitude only. 
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4.3 Set of Variational Equations for Numerical 
Investigations 

4.3.1 UNCONSTRAINED ARCS 


From equations (4.1) to (4.9) follows, with the Hamiltonian function 


H = - A g sin y + — cos y A (1 - gr) - X cos x tan A 
v r ' [ 7 X 


COS Y 

+ X. — + X. sin x I + A, A sin y 

A cos A A | h 


+ &L V C-' Sh 
2m 


M 


illl xl 


X SOUL + 1 COB M 

X cos y A 


) 


VX c „ 
v W 


, (4.3.1) 


a set of differential equations for the Lagrangian multipliers according to 
X = -H 


x 


0 


cos y 


cos X 
cos A 


(\ slnA - A e) - \ ( 1+ v£) 


- A a sin X 


2 VA c - c. 


. ^ S -/3h 

A, sin y + r-“* e 
h 2m 


. _ - . sin M + A cos M 

v D L\cosy y , 


A = 


V cos y 
r 


c 


e 


\ sin x 

A sinA - A, cos X 

X / cos A A 


(4.3.2) 


(4.3.3) 
, (4.3.4) 


A = 


A g 0 R 
cos y \ — 


A, V ) + 
n 


x 


cos X /. . . 

I A “ A sin 

cos A \ 0 X 


111 A ) + \ \ 1 - 


+ A a sin X 


A sin y ~ 

— — ~ ^ e ^ c V sin P 
cosV 2m L 


(4.3.5) 
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(4.3.6) 


A. = 


A, = 


V cos y cos x / -v . A \ 

S [ A - A sinAl , 

r cos A y X 9 J 

2 g n R 2 l . y i V cos y 

- yA^ sin y + -f cos yj + 

CQS - (a - A sin A ) + A + A. sin x 
cos A \ 9 X ) y A 


x 


VA c^ - c T I sin P + A cos P 

v D L \ cos y y 


2m 


(4.3.7) 


This set is valid for unconstrained arcs with the controls given by 


p = tan 
opt 


-1 


A /[A cosy\ 

V V r y 


(4.3.8) 


and 


-w 


l/(n-l) 


(4.3.9) 


L opt lVA v kn 


where 



4.3.2 CONSTRAINED ARCS 


On a constrained arc of the form 

c = h - h e ( C L , v ) = 0 , 

Reference 25, the control is chosen by satisfying the constraint 


(4.3.10) 


(4.3.11) 


and there appear additive terms in the adjoint differential equations (second 
summand) 
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(4.3.12) 


Since the only control involved in C is the lift coefficient C and the 

-Li 

-ates involved are h and V , one obtains 


8C 


X Y A v + K 9V 

c uc 


LH 


(4.3.13) 


0C 


A., — A, - K 

h h, uc 3h 

c 


LH 


(4.3.14) 


with 


K = 


spy 

2m 


\ (n-l) 

sin P + A cos P - V A kn C T „ 

cos y y v LH 


C and its partials are determined according to the model described in 
LH 

section 3. 4. 


4.3.3 BOUNDARY CONDITIONS 

The boundary conditions to these and the state differential equations, 
equations (4. 1. 1) to (4. 1. 9) , are as follows: The initial state variables 
are considered to be given. At the final time the altitude is fixed, h^ = 30 1cm 

was taken as stopping condition for the integration. With all other state 
variables open to assume the most favorable values for maximum range, the 
trajectories end with a flare, trading kinetic energy for range in the final 
part and the lift coefficient goes to its maximum value [ 24] . This undesirable 
boundary layer type of behavior was suppressed by prescribing vertical 
equilibrium at the final time with a lift coefficient corresponding to maximum 
L/D . For open final heading /a = o\ it is seen from equation (4. 3. 8) 

\ *f / 

that the bank angle is P^ = 0 . With altitude, lift coefficient, 
given, equation (4. 1. 3) can be solved for the QSG velocity 
using equation (4. 1.42a) to eliminate the flight path angle y ^ . 


and bank angle 

(r f = o) 
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However, to reduce the number of adjoint sets to be integrated back- 
ward in the min-H-gradient procedure [ 24] , and y^ were not prescribed 

directly. In a quasisteady glide state, which is desirable for continuation of 
the trajectory to landing approach, the Lagrangian multiplier to the velocity 
is easily determinable from the velocity by a physical consideration. The 
aerodynamic properties of the vehicle translate the kinetic energy altitude 



into range via the factor E = L/D . Therefore, for range R measured in 
a geocentric angle, there follows for the conditions given at final time 



(4.3.16) 


With V and 
equation (4.3. 


A given, the multiplier A can now be determined from 

V f T f 

9) for the desired final lift coefficient C : 

L.E 




kn C 


(n-1) 


LE 


(4.3.17) 


For general range maximization, the payoff quantity is 

<fi = A + p© . (4.3.18) 

p = 0 leads to maximum lateral range A and p » 1 leads to maximum 

UlcLX 

downrange 0 . Therefore, one can expect differences between optimal 

max 

s 

trajectories only for values of p in orders of magnitude around 1 (p = 10 
with j s | < say 3 ) [ 24] . This payoff quantity leads to final values for the 
Lagrangian multipliers. 

= -1 


= A Q = -p . (4.3.19) 
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For footprint evaluations with QSG final conditions, the boundary value prob- 
lem has, thus, been reduced to requiring only the backward integration of 
one set of adjoint equations without resorting to penally functions. 


5. HEATING CONSTRAINED MAXIMUM LATERAL 
RANGE NON-QSG ENTRY TRAJECTORY 


For the heating constraint modeled in section 3. 4 with AC = 0. 03 , 

a maximum lateral range trajectory for a Space Shuttle orbiter-type vehicle 
has been determined (Fig. 7). 

The numerical iteration with a refined gradient program was stopped 
when the increase in lateral range had become very slow. Of the control 
update computed from the min-H -feature, only about 1 percent could be used 
without having the iteration diverge. Figure 7d shows altitude versus velocity 
for three cases: The dash-dotted smooth curve is a QSG trajectory for 
comparison, around which the (solid) unconstrained trajectory with V =7.850 

km/s at 80 km altitude and y = -1.25 deg is seen to oscillate. There are 

G 

five peaks in the velocity range down to 4 km/s. The heating constrained 
trajectory has only three peaks in this region and is shifted upward. This is 
achieved by a higher lift coefficient (Fig. 7c) and a smaller bank angle (Fig. 

7b) , initially. When the heating constraint is first encountered, the lift 
coefficient yields in order not to violate the constraint. The iteration process 
positioned this dip around the lift coefficient for best L/D . During the 
first dip , the bank angle is slightly increased to use the lift for a heading 
change and decrease the upward velocity which would lead to too large a flight 
path angle at the beginning of the second constrained arc where the heating 
limit is most severe due to the combined effect of high velocity and beginning 
gravity pull, no more balanced by the centrifugal force. In this region, the 
bank angle decreases below the value for the unconstrained trajectory and 
stays there for a while. After the heating region is passed, the bank angle 
is increased above that for the unconstrained trajectory to make up for the 
heading change loss in the earlier parts. As Figure 7a shows, this results 
in a downrange increase of approximately 4 deg, but the lateral range loss is 
only about 2 percent. The arcs on the heating boundary are approximately 
4^ , 7, and 2 5 minutes long. 
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6. CONCLUSIONS 


Maximum range, three-dimensional atmospheric entry trajectories 
for lifting vehicles have been approximated by a superposition of a non- 
oscillatory quasisteady glide component and solutions to linearized perturba- 
tion equations. For the QSG motion, variational equations have been derived. 
They show that the multiplier to the flight path angle is a linear combination 
of the other multipliers, which may be the cause of numerical difficulties 
encountered when applying indirect methods to range maximization. 

By applying the maximum principle to the QSG problem, a simple 
bank angle control law has been obtained without taking the boundary conditions 
into account. This "basic” control, in comparing it to numerical solutions 
of the full equations, turns out to be closely followed for maximum lateral 
range trajectories. Adaptation to the special case is achieved by modifications 
propagating from the final time backward. Because of the simplicity of this 
basic control, bank angle equal to arc sin of ratio velocity to local satellite 
velocity and nondimensional bank angle rate approximately equal to inverse of 
lift-to-drag ratio, some analytical approximate solutions have been obtained. 
The maximum cross-range maneuver time can be determined fairly accurately. 
A full analytical solution, however was not achieved. Numerical QSG solu- 
tions are given for Space Shuttle orbiter class vehicles with L/D from 1.4 
to 3. For general footprint trajectories, the optimal bank angle program is 
not linear with time. The perturbation equations around QSG show that the 
frozen state oscillation period is almost independent of vehicle parameters 
and depends only on velocity while the damping decreases with increasing 
L/D . There is an unstable motion component over the entire velocity range, 
increasing rapidly as v becomes less than 1/3. In this region the damping 
and the frequency of the oscillatory mode also increase strongly, leading to the 
almost smooth altitude-velocity curves for all perturbed trajectories at these 
low velocities. The range of validity of the simple linear perturbation equa- 
tions is limited to a few tenths of a degree in flight path angle and a few 
kilometers in altitude. 

The influence of an angle of attack (lift coefficient) , altitude- and 
velocity-dependent heating constraint on a maximum lateral range earth- 
orbit entry trajectory for a Space Shuttle orbiter-type vehicle has been inves- 
tigated numerically. Although the shape of the control time histories and the 
trajectory have been changed considerably, the loss in lateral range was 
only about 2 percent compared to the unconstrained one. 

George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, May 1972 
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